function schur_decomposition()
% SCHUR_DECOMPOSITION Schur分解演示
% 
% 展示Schur分解的理论、计算和应用

fprintf('=== Schur分解演示 ===\n\n');

%% 1. Schur分解基础
fprintf('1. Schur分解基础\n');
fprintf('目标：理解Schur分解的基本概念\n\n');

% 简单的3x3矩阵
A1 = [4, 1, 2; 0, 3, 1; 0, 0, 2];
fprintf('上三角矩阵 A1 =\n');
disp(A1);

% 计算Schur分解
[U1, T1] = schur(A1);
fprintf('正交矩阵 U1 =\n');
disp(U1);
fprintf('Schur形式 T1 =\n');
disp(T1);

% 验证分解
A1_reconstructed = U1 * T1 * U1';
fprintf('重构误差: %.2e\n', norm(A1 - A1_reconstructed, 'fro'));

% 验证U的正交性
orthogonality1 = norm(U1' * U1 - eye(size(U1)), 'fro');
fprintf('U的正交性误差: %.2e\n', orthogonality1);

% 特征值在对角线上
eigenvals1 = diag(T1);
eigenvals1_eig = eig(A1);
fprintf('Schur形式对角元: [%.4f, %.4f, %.4f]\n', eigenvals1);
fprintf('eig函数结果: [%.4f, %.4f, %.4f]\n', sort(eigenvals1_eig, 'descend'));

%% 2. 一般矩阵的Schur分解
fprintf('\n2. 一般矩阵的Schur分解\n');
fprintf('目标：处理一般非对称矩阵\n\n');

A2 = [1, 2, 3; 4, 5, 6; 7, 8, 9];
fprintf('一般矩阵 A2 =\n');
disp(A2);

[U2, T2] = schur(A2);
fprintf('Schur形式 T2 =\n');
disp(T2);

% 检查上三角结构
lower_part = tril(T2, -1);
fprintf('下三角部分的范数: %.2e (应该接近0)\n', norm(lower_part, 'fro'));

% 特征值
eigenvals2 = diag(T2);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals2);

% 验证重构
A2_reconstructed = U2 * T2 * U2';
fprintf('重构误差: %.2e\n', norm(A2 - A2_reconstructed, 'fro'));

%% 3. 复Schur分解 vs 实Schur分解
fprintf('\n3. 复Schur分解 vs 实Schur分解\n');
fprintf('目标：比较复数和实数Schur分解\n\n');

% 有复特征值的矩阵
A3 = [0, 1; -1, 0];  % 旋转矩阵
fprintf('旋转矩阵 A3 =\n');
disp(A3);

% 复Schur分解
[U3_complex, T3_complex] = schur(A3, 'complex');
fprintf('复Schur形式 T3_complex =\n');
disp(T3_complex);

% 实Schur分解
[U3_real, T3_real] = schur(A3, 'real');
fprintf('实Schur形式 T3_real =\n');
disp(T3_real);

% 分析实Schur形式的2x2块
fprintf('实Schur形式包含2×2块，对应复共轭特征值对\n');
eigenvals3 = eig(A3);
fprintf('特征值: %.4f ± %.4fi\n', real(eigenvals3(1)), abs(imag(eigenvals3(1))));

% 从2x2块提取特征值
block_2x2 = T3_real;
eigenvals_from_block = eig(block_2x2);
fprintf('从2×2块计算的特征值: %.4f ± %.4fi\n', ...
        real(eigenvals_from_block(1)), abs(imag(eigenvals_from_block(1))));

%% 4. 有序Schur分解
fprintf('\n4. 有序Schur分解\n');
fprintf('目标：按特征值大小排序Schur分解\n\n');

A4 = [3, 1, 4; 1, 5, 9; 2, 6, 5];
fprintf('矩阵 A4 =\n');
disp(A4);

% 标准Schur分解
[U4, T4] = schur(A4);
eigenvals4_unsorted = diag(T4);
fprintf('未排序的特征值: [%.4f, %.4f, %.4f]\n', eigenvals4_unsorted);

% 有序Schur分解（按模长降序）
[U4_ordered, T4_ordered] = ordschur(U4, T4, 'lhp');  % 左半平面优先
eigenvals4_ordered = diag(T4_ordered);
fprintf('排序后的特征值: [%.4f, %.4f, %.4f]\n', eigenvals4_ordered);

% 验证排序
[~, sort_idx] = sort(abs(eigenvals4_unsorted), 'descend');
fprintf('预期排序: [%.4f, %.4f, %.4f]\n', eigenvals4_unsorted(sort_idx));

%% 5. Schur分解的应用：矩阵函数计算
fprintf('\n5. Schur分解的应用：矩阵函数计算\n');
fprintf('目标：使用Schur分解计算矩阵函数\n\n');

A5 = [2, 1; 0, 3];
fprintf('矩阵 A5 =\n');
disp(A5);

[U5, T5] = schur(A5);

% 计算矩阵指数 exp(A)
% 对于上三角矩阵，可以直接计算
T5_exp = zeros(size(T5));
n = size(T5, 1);

% 对角元素
for i = 1:n
    T5_exp(i, i) = exp(T5(i, i));
end

% 上三角元素（使用递推公式）
for i = 1:n-1
    for j = i+1:n
        if abs(T5(i,i) - T5(j,j)) > 1e-12
            T5_exp(i, j) = T5(i, j) * (T5_exp(j, j) - T5_exp(i, i)) / (T5(j, j) - T5(i, i));
        else
            T5_exp(i, j) = T5(i, j) * T5_exp(i, i);  % 重特征值情况
        end
    end
end

A5_exp = U5 * T5_exp * U5';
fprintf('通过Schur分解计算的 exp(A5) =\n');
disp(A5_exp);

% 与MATLAB内置函数比较
A5_exp_matlab = expm(A5);
fprintf('MATLAB expm结果 =\n');
disp(A5_exp_matlab);
fprintf('差异: %.2e\n', norm(A5_exp - A5_exp_matlab, 'fro'));

%% 6. 广义Schur分解
fprintf('\n6. 广义Schur分解\n');
fprintf('目标：求解广义特征值问题的Schur形式\n\n');

A6 = [1, 2; 3, 4];
B6 = [2, 1; 1, 3];
fprintf('矩阵 A6 =\n'); disp(A6);
fprintf('矩阵 B6 =\n'); disp(B6);

% 广义Schur分解：A = Q*S*Z', B = Q*T*Z'
[AA, BB, Q6, Z6] = qz(A6, B6);
fprintf('广义Schur形式:\n');
fprintf('S (对应A) =\n'); disp(AA);
fprintf('T (对应B) =\n'); disp(BB);

% 广义特征值 = diag(S) ./ diag(T)
generalized_eigenvals = diag(AA) ./ diag(BB);
fprintf('广义特征值: [%.6f, %.6f]\n', generalized_eigenvals);

% 验证
eigenvals6_eig = eig(A6, B6);
fprintf('eig(A,B)结果: [%.6f, %.6f]\n', eigenvals6_eig);
fprintf('差异: %.2e\n', norm(sort(generalized_eigenvals) - sort(eigenvals6_eig)));

%% 7. Schur分解的数值稳定性
fprintf('\n7. Schur分解的数值稳定性\n');
fprintf('目标：分析Schur分解的数值特性\n\n');

% 测试不同条件数的矩阵
fprintf('条件数        重构误差      正交性误差\n');
fprintf('--------      --------      ----------\n');

condition_numbers = [1e2, 1e6, 1e10, 1e14];
for cond_num = condition_numbers
    % 构造指定条件数的矩阵
    n = 4;
    eigenvals = [1, 1/sqrt(cond_num), 1/cond_num^(1/3), 1/cond_num];
    V_rand = orth(randn(n, n));
    A_test = V_rand * diag(eigenvals) * V_rand';
    
    % Schur分解
    [U_test, T_test] = schur(A_test);
    
    % 评估误差
    reconstruction_error = norm(A_test - U_test * T_test * U_test', 'fro') / norm(A_test, 'fro');
    orthogonality_error = norm(U_test' * U_test - eye(n), 'fro');
    
    fprintf('%.0e      %.2e        %.2e\n', cond_num, reconstruction_error, orthogonality_error);
end

%% 8. Schur分解 vs 特征分解
fprintf('\n8. Schur分解 vs 特征分解\n');
fprintf('目标：比较两种分解的特点\n\n');

% 非对称矩阵
A8 = [1, 2, 3; 0, 4, 5; 0, 0, 6];
fprintf('上三角矩阵 A8 =\n');
disp(A8);

% Schur分解
tic;
[U8, T8] = schur(A8);
time_schur = toc;

% 特征分解
tic;
[V8, D8] = eig(A8);
time_eigen = toc;

fprintf('计算时间比较:\n');
fprintf('Schur分解: %.6f 秒\n', time_schur);
fprintf('特征分解: %.6f 秒\n', time_eigen);

% 数值稳定性比较
cond_U = cond(U8);
cond_V = cond(V8);
fprintf('条件数比较:\n');
fprintf('Schur分解的U矩阵: %.2e\n', cond_U);
fprintf('特征分解的V矩阵: %.2e\n', cond_V);

fprintf('Schur分解的优势:\n');
fprintf('• U矩阵总是正交的（条件数=1）\n');
fprintf('• 数值稳定性更好\n');
fprintf('• 适用于所有矩阵\n');
fprintf('• 便于计算矩阵函数\n\n');

%% 总结
fprintf('=== Schur分解总结 ===\n');
fprintf('1. Schur分解将矩阵化为上三角形式\n');
fprintf('2. 正交变换矩阵保证数值稳定性\n');
fprintf('3. 复Schur分解得到复上三角矩阵\n');
fprintf('4. 实Schur分解保持实数运算\n');
fprintf('5. 可以按特征值排序\n');
fprintf('6. 适用于矩阵函数计算\n');
fprintf('7. 广义Schur分解处理矩阵束\n');
fprintf('8. 比特征分解更稳定\n\n');

end